{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 灰色预测模型(Grey Model, GM)\n",
    "\n",
    "**只适合中短期预测**\n",
    "\n",
    "1. GM(1, 1)建模\n",
    "    - 要求: 需要一次累加序列近似为指数规律发展\n",
    "    - 模型使用数据不是原有数据\n",
    "\n",
    "2. GM(2, 1), DGM, Verhulst\n",
    "    - 对于非单调的摆动发展序列或饱和的S形序列.\n",
    "\n",
    "3. 灰色Verhulst预测模型\n",
    "    - 用来描述具有饱和状态的S型模型. 常用于人口预测, 生物生长, 繁殖预测, 产品经济寿命预测等.\n",
    "    \n",
    "- 已知$x^{(0)}=(x^{(0)}(1),x^{(0)}(2),\\cdots,x^{(0)}(n))$\n",
    "- 累加生成序列:$x^{(1)}=(x^{(0)}(1),x^{(0)}(1)+x^{(0)}(2),\\cdots,x^{(0)}(1)+\\cdots+x^{(0)}(n))$\n",
    "- 均值生成序列:$z^{(1)}=(z^{(1)}(1),z^{(1)}(2),\\cdots,z^{(1)}(n))$, 其中$z^{(1)}(k)=0.5x^{(1)}(k)+0.5x^{(1)}(k-1)$\n",
    "- 灰微分方程: $x^{(0)}(k)+az^{(1)}(k)=b,\\ k=2,3,\\cdots,n$\n",
    "- 相应白化微分方程: $\\frac{dx^{(1)}}{dt}+ax^{(1)}(t)=b$\n",
    "- 可容覆盖: $\\Theta=(e^{-\\frac{2}{n+1}}, e^{\\frac{2}{n+1}})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "级比: [0.9820442  1.         1.00416089 1.00980392 0.99166667 1.00558659]\n",
      "可容覆盖: (0.7788007830714049, 1.2488488690016821)\n",
      "可容性: [True, True, True, True, True, True]\n",
      "一次累加序列: [ 71.1 143.5 215.9 288.  359.4 431.4 503. ]\n",
      "B矩阵:\n",
      "[[-107.3     1.  ]\n",
      " [-179.7     1.  ]\n",
      " [-251.95    1.  ]\n",
      " [-323.7     1.  ]\n",
      " [-395.4     1.  ]\n",
      " [-467.2     1.  ]]\n",
      "Y值: [72.4 72.4 72.1 71.4 72.  71.6]\n",
      "最小二乘法得到a, b: [2.34378648e-03 7.26572696e+01]\n",
      "x1预测值: [ 71.1        143.50574144 215.741978   287.8091065  359.70752283\n",
      " 431.43762195 502.9997979 ]\n",
      "x0预测值: [71.1        72.40574144 72.23623656 72.0671285  71.89841633 71.73009912\n",
      " 71.56217595]\n",
      "\n",
      "   Year Raw Data Predction    Error Relative_Error     级比偏差\n",
      "0  1986     71.1   71.1000   0.0000          0.00%      nan\n",
      "1  1987     72.4   72.4057  -0.0057         -0.01%   0.0203\n",
      "2  1988     72.4   72.2362   0.1638          0.23%   0.0023\n",
      "3  1989     72.1   72.0671   0.0329          0.05%  -0.0018\n",
      "4  1990     71.4   71.8984  -0.4984         -0.70%  -0.0074\n",
      "5  1991     72.0   71.7301   0.2699          0.37%   0.0107\n",
      "6  1992     71.6   71.5622   0.0378          0.05%  -0.0032\n"
     ]
    }
   ],
   "source": [
    "# GM(1, 1)建模\n",
    "# 要求: 需要一次累加序列近似为指数规律发展\n",
    "# 例15.2\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# 1.级比检验\n",
    "x0 = np.array([71.1, 72.4, 72.4, 72.1, 71.4, 72.0, 71.6], dtype=np.float)\n",
    "lamb = np.array([x0[i] / x0[i+1] for i in range(x0.shape[0]-1)])\n",
    "print('级比:', lamb)\n",
    "Theta = (np.exp(-2/(x0.shape[0]+1)), np.exp(2/(x0.shape[0]+2)))\n",
    "print('可容覆盖:', Theta)\n",
    "OK = [i <= Theta[1] and i >= -Theta[0] for i in lamb]\n",
    "print('可容性:', OK)\n",
    "\n",
    "# GM(1, 1)建模\n",
    "x1 = np.array([x0[:i].sum() for i in range(1, x0.shape[0]+1)])\n",
    "print('一次累加序列:', x1)\n",
    "B = np.array([[-(x1[i]+x1[i+1])/2, 1] for i in range(x1.shape[0]-1)])\n",
    "Y = np.array([x0[i] for i in range(1, x0.shape[0])])\n",
    "print('B矩阵:', B, sep='\\n')\n",
    "print('Y值:', Y)\n",
    "u = np.matmul(np.linalg.inv(np.matmul(B.T, B)), np.matmul(B.T, Y))\n",
    "print('最小二乘法得到a, b:', u)\n",
    "a, b = u[0], u[1]\n",
    "x1_k = lambda k: (x0[0] - b/a) * np.exp(-a*k) + b/a\n",
    "\n",
    "# 预测\n",
    "x1_pred = np.array([x1_k(k) for k in range(x0.shape[0])])\n",
    "print('x1预测值:', x1_pred)\n",
    "x0_pred = np.array([x1_pred[i] - x1_pred[i-1] if i != 0 else x1_pred[0] for i in range(x1_pred.shape[0])])\n",
    "print('x0预测值:', x0_pred)\n",
    "\n",
    "# 检验\n",
    "test_table_columns = ['Year', 'Raw Data', 'Predction', 'Error', 'Relative_Error', '级比偏差']\n",
    "test_table = pd.DataFrame(np.c_[[str(i) for i in range(1986, 1993)],\n",
    "                                x0,\n",
    "                                x0_pred,\n",
    "                                x0 - x0_pred,\n",
    "                                np.divide((x0 - x0_pred), x0),\n",
    "                                [np.nan] + [1-(1-0.5*a)/(1+0.5*a)*lamb[k] for k in range(lamb.shape[0])]],\n",
    "                          columns=test_table_columns)\n",
    "test_table['Predction'] = test_table['Predction'].apply(lambda x: format(float(x), '.4f'))\n",
    "test_table['Error'] = test_table['Error'].apply(lambda x: format(float(x), '.4f'))\n",
    "test_table['Relative_Error'] = test_table['Relative_Error'].apply(lambda x: format(float(x), '.2%'))\n",
    "test_table['级比偏差'] = test_table['级比偏差'].apply(lambda x: format(float(x), '.4f'))\n",
    "print('\\n', test_table, sep='')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "级比: [0.83673469 0.80327869 0.78205128 0.8125     0.92307692]\n",
      "可容覆盖: (0.751477293075286, 1.2840254166877414)\n",
      "可容性: [True, True, True, True, True]\n",
      "1-AGO序列: [ 41.  90. 151. 229. 325. 429.]\n",
      "1-IAGO序列: [ 8. 12. 17. 18.  8.]\n",
      "均值生成序列: [ 65.5 120.5 190.  277.  377. ]\n",
      "B矩阵:\n",
      "[[ -49.   -65.5    1. ]\n",
      " [ -61.  -120.5    1. ]\n",
      " [ -78.  -190.     1. ]\n",
      " [ -96.  -277.     1. ]\n",
      " [-104.  -377.     1. ]]\n",
      "Y值: [ 8. 12. 17. 18.  8.]\n",
      "最小二乘法得到a1, a2, b: [ -1.09219635   0.19590335 -31.79834712]\n",
      "得到带参符号解: Eq(x1_t(t), C1*exp(t*(-a1 - sqrt(a1**2 - 4*a2))/2) + C2*exp(t*(-a1 + sqrt(a1**2 - 4*a2))/2) + b/a2)\n",
      "待定常数:  {C1, C2}\n",
      "关于待定常数的多元方程组:\n",
      "[Eq(41, C1 + C2 + b/a2), Eq(429.0, C1*exp(-5*a1/2 - 5*sqrt(a1**2 - 4*a2)/2) + C2*exp(-5*a1/2 + 5*sqrt(a1**2 - 4*a2)/2) + b/a2)]\n",
      "解得待定常数: {C1: ((41.0*a2 - b)*exp(2.5*a1 + 5.0*sqrt(a1**2 - 4.0*a2)) - (429.0*a2 - b)*exp(5.0*a1 + 2.5*sqrt(a1**2 - 4.0*a2)))*exp(-2.5*a1)/(a2*(exp(5.0*sqrt(a1**2 - 4.0*a2)) - 1.0)), C2: -((41.0*a2 - b)*exp(2.5*a1) - (429.0*a2 - b)*exp(5.0*a1 + 2.5*sqrt(a1**2 - 4.0*a2)))*exp(-2.5*a1)/(a2*(exp(5.0*sqrt(a1**2 - 4.0*a2)) - 1.0))} \n",
      "\n",
      "最终解: x1_t(t) = Eq(x1_t(t), 203.84901286639*exp(0.226223404169046*t) - 0.532505769427929*exp(0.865972945416781*t) - 162.316507096962) \n",
      "\n",
      "x1预测值: [41.0000000000000 92.0148144078681 155.156052197627 232.367192369488\n",
      " 324.521982077495 429.000000000000]\n",
      "x0预测值: [41.0000000000000 51.0148144078680 63.1412377897586 77.2111401718614\n",
      " 92.1547897080069 104.478017922505]\n",
      "\n",
      "  Raw Data Predction    Error Relative_Error\n",
      "0       41   41.0000  -0.0000         -0.00%\n",
      "1       49   51.0148  -2.0148         -4.11%\n",
      "2       61   63.1412  -2.1412         -3.51%\n",
      "3       78   77.2111   0.7889          1.01%\n",
      "4       96   92.1548   3.8452          4.01%\n",
      "5      104  104.4780  -0.4780         -0.46%\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZzNdRfA8c9pDGPfwuOhB0URjcEknpCKUtbqmZCkUlJ2sraXkFAkSaJJsmTJvjW2QZaZUNZQZDIxlH2Z7fv88b00psGYe2d+994579fL68793d+de6Zy+s75nd/5ijEGpZRS/uUGpwNQSinleZrclVLKD2lyV0opP6TJXSml/JAmd6WU8kM5nA4A4MYbbzRly5Z1OgyllPIp0dHRR40xxdJ6zSuSe9myZYmKinI6DKWU8ikicuBKr2lZRiml/JAmd6WU8kOa3JVSyg95Rc09LQkJCcTExHD+/HmnQ/EbQUFBlC5dmsDAQKdDUUplMq9N7jExMeTPn5+yZcsiIk6H4/OMMRw7doyYmBjKlSvndDhKqUzmtWWZ8+fPU7RoUU3sHiIiFC1aVH8TUiqb8NrkDmhi9zD956lU9uG1ZRmllHJXcjIkJf39JzHx8udXO57eY+6+v3JlePxxz//smtyvIiAggDvuuIPExETKlSvHpEmTKFSokMe+f/369YmNjSVXrlzEx8fToEEDBg4ceM3PGDRoEAMGDPBYHEr5ivh4ePVVmD4dEhKunXB9QcuWmtyzXO7cudmyZQsA7dq14+OPP+aVV17x6GdMnjyZ0NBQ4uPj6d+/P82bN2fVqlVXfY8md5UdHThgE+GGDdAs5DeK1fgPOXJAQMDlf9I6dqXjTp97ww2QWdVSTe7pVLt2bX788UcATp8+TfPmzfnrr79ISEhg4MCBNG/enKFDhxIUFETXrl3p0aMHW7duZfny5URERDBx4kS++uqrK37/nDlzMnToUMqXL8/WrVupWrUqLVq04ODBg5w/f55u3brRoUMH+vXrx7lz5wgJCaFy5cpMnjw5zfOU8ifz58NTbRJJOhfPN7TjfwciYM1ByJvX6dC8lm8k9+7dwbWC9piQEPjww3SdmpSUREREBO3btwdsv/js2bMpUKAAR48epVatWjRr1ox69eoxfPhwunbtSlRUFBcuXCAhIYE1a9ZQt27da35OQEAAVatWZdeuXVStWpUJEyZQpEgRzp07x5133sljjz3GkCFDGD169KXfKIA0zytatGjG/rko5UUSEuCV5w7z/pclqMaPfJPnaW55uQn0GKOJ/Rp8I7k75OIKef/+/dSoUYOGDRsCtmd8wIABrF69mhtuuIHff/+dw4cPU6NGDaKjozl16hS5cuWievXqREVFERkZyahRo9L1mSn3tB01ahSzZ88G4ODBg+zZsyfNpJ3e85TyJQfnbaFVu5ys++t2Xsw5nhEvxxLUayUUKeJ0aD7BN5J7OlfYnnax5n7ixAmaNGnCxx9/TNeuXZk8eTJxcXFER0cTGBhI2bJlOX/+/KWvJ06cyH//+1+Cg4NZsWIF+/bto1KlStf8vKSkJH766ScqVarEypUr+e677/j+++/JkycP9evXT7NHPb3nKeUz1q5lUddFtP2hOxcIYkrYLFp9FgYFCzodmU/x6j53b1GwYEFGjRrFsGHDSEhI4MSJExQvXpzAwEBWrFjBgQN/T92sV68ew4YNo169etStW5exY8cSEhJyzR7zhIQE+vfvz0033URwcDAnTpygcOHC5MmTh127drF+/fpL5wYGBpKQkABw1fOU8hnGwMqVJNZvwIA6q3j4h4GUKplMdDS0mv6oJvYM0OSeTtWqVaNq1apMnTqVNm3aEBUVRWhoKJMnT6ZixYqXzqtbty6xsbHUrl2bEiVKEBQUdNV6e5s2bQgODqZKlSqcOXOGOXPmANCoUSMSExMJDg7mtddeo1atWpfe06FDB4KDg2nTps1Vz1PK6xkDS5dCvXocuvcJ7l/3DoMZwPNPJ7B+X3FurZ7P6Qh9lqSs8TolNDTUpN6sY+fOnekqZajro/9clVcwBhYuhHfegQ0bWHZja9qcH8+Z5Nx8+qnw5JNOB+gbRCTaGBOa1mu6cldKZZ3kZJg9G0JDoUkTkv6I4/XGUTx4bDLFy+QhKkoTu6docldKZb6kJHtbaUgIPPoonDzJHx9MoWG5PbyzoAZPPy1s3Aj6S6Xn+Ea3jFLKNyUmwrRpMHAg7NoFFSvCV1+xvHgrnmgbwMmTMHEiPP2004H6n2uu3EVkgogcEZFtKY4VEZFlIrLH9Vg4xWv9RWSviOwWkQczK3CllBdLSLBZu1IlePJJCAyEadNI2rqNt/e1oWGjAAoXho0bNbFnlvSUZb4AGqU61g+IMMZUACJczxGR24FWQGXXe8aISIDHolVKebcLF2DcOLj1Vnj2WcifH2bNgi1bOFL/cRo1DuCNN+CJJ2DTJqhSxemA/dc1k7sxZjXwZ6rDzYFw19fhQIsUx6caYy4YY34F9gI1PRSrUspbnT8Po0dD+fLwwgtQvLgdCBMdDY88wqrIGwgJgTVr4LPP4MsvIZ92OWaqjF5QLWGMiQVwPRZ3HS8FHExxXozr2D+ISAcRiRKRqLi4uAyGkbkCAgIICQmhSpUqhIWFcfbs2Qx/r5UrV9KkSZM0jxcsWJBq1apx2223Ua9ePebPn5+u77du3boMx6OUR5w9Cx98ADffDF26QNmysGQJrF8PjRuTbIRBg+C++2wy37ABnnsu8yYhqr95ulsmrX9laTbSG2PGGWNCjTGhxYoV83AYnnFx/MC2bdvImTMnY8eOvex1YwzJycluf07dunXZvHkzu3fvZtSoUXTu3JmIiIirvkeTu3LU6dMwdCiUKwc9e9oLpStWwOrV8MADIMLRo9C4Mbzyip1XHh0NwcFOB559ZDS5HxaRkgCuxyOu4zHATSnOKw0cynh43qNu3brs3buX/fv3U6lSJV566SWqV6/OwYMHWbp0KbVr16Z69eqEhYVx+vRpABYvXkzFihWpU6cOs2bNStfnhISE8PrrrzN69GgA5s2bx1133UW1atVo0KABhw8fZv/+/YwdO5YPPviAkJAQIiMj0zxPKY87cQLefRfKlIG+faFaNVtrWb4c6te/tCRfu9Z2PS5fDp98Al9/bcvvKutktBVyLtAOGOJ6nJPi+NciMgL4N1AB2OhukA5P/CUxMZFFixbRqJG9rrx7924mTpzImDFjOHr0KAMHDuS7774jb968vPfee4wYMYI+ffrw/PPPs3z5csqXL0/Lli3THVv16tV5//33AahTpw7r169HRBg/fjxDhw5l+PDhdOzYkXz58vHyyy8D8Ndff6V5nlIe8eefMHIkjBoFx49Dkybw2mtQ8/JLasnJMGwYDBhgKzTr19v8r7LeNZO7iEwB6gM3ikgM8AY2qU8XkfbAb0AYgDFmu4hMB3YAiUAnY0xSJsWe6S6O/AW7cm/fvj2HDh2iTJkyl2a4rF+/nh07dnD33XcDEB8fT+3atdm1axflypWjQoUKADz55JOMGzcuXZ+bciRETEwMLVu2JDY2lvj4eMqVK5fme9J7nlLX5ehRGDHCXiw9dQoeecTuc1e9+j9OPXYM2rWDBQvgf/+D8eN13peTrpncjTGtr/DS/Vc4/13gXXeCSs2hib+XbbOXUt4UmwQYY2jYsCFTpky57JwtW7ZccxLklWzevPnS/JcuXbrQs2dPmjVrxsqVK3nzzTfTfE96z1MqXf74A4YPhzFj4Nw5WzR/5RW44440T1+/3m6BFxsLH30EnTrpRVOn6fgBN9WqVYu1a9eyd+9eAM6ePcvPP/9MxYoV+fXXX9m3bx/AP5L/lfz444+88847dOrUCbAjfUuVsg1H4eHhl87Lnz8/p06duvT8SucpdV1+/x26dbMXSkeMsKMCtm+HqVPTTOzG2GaZunXtfqDr1kHnzprYvYEmdzcVK1aML774gtatWxMcHEytWrXYtWsXQUFBjBs3jsaNG1OnTh3KlClzxe8RGRl5qRWyU6dOjBo1ivvvt78Yvfnmm4SFhVG3bl1uvPHGS+9p2rQps2fPvnRB9UrnKZUuBw7ASy/ZlsYxY+xdRrt3w6RJVxz48tdftkrTs6ftivnhBzsPTHkHHfmbzeg/V3WZfftg8GAID7fL7WefhX797NXQq9i0yVZqYmJsR2T37rpad8LVRv7q4DClsqPdu2HQIJg8GXLkgI4doU8fuOmmq77NGHtttVcvKFkSIiNB94fxTprclcpOtm+3ferTpkGuXLa+/vLLNlNfw4kT0L49zJxpOyHDw3Wvam/m1cndGJPhjhP1T95QglMOiYuzNfUZM+wcgN69bbG8ePFrvxdbT3/8cdi/35ZhevWyF1CV9/La5B4UFMSxY8coWrSoJngPMMZw7NgxgoKCnA5FZbVTp+Chh+yq/dVXbYG8aNF0vdUYGDvWvqVYMVi1Cly3dCgv57XJvXTp0sTExOCtQ8V8UVBQEKVLl3Y6DJWV4uNtO+OWLTBnjm1rSaeTJ6FDB1vBadTINs5oI5bv8NrkHhgYqHdZKuWO5GR7y+h338EXX1xXYt+6FcLCbDPNoEF2jIyWYXyL1yZ3pZQbjIEePezNR0OH2iSfzreNHw9du0LhwnbQY716mRyryhT6/2Kl/NGQIXbIV8+ethsmHU6fhrZtbSmmTh1bydHE7rs0uSvlbz7/3I5lfPJJeP/9dN1dtG0b3HmnHc371luweHG6G2mUl9KyjFL+ZO5cu/Ru1AgmTEhXoXziRDvoq0ABW56/774siFNlOl25K+UvIiPtaMbQUPjmGwgMvOrpZ87A00/biQO1atkyjCZ2/6HJXSl/8NNP0KyZ3SFpwYJr7j69c6fdZ+PLL+2eG8uWwb/+lUWxqiyhZRmlfN2BA7YMkyeP3Zz6Gs3okybZUTJ589ra+gMPZFGcKkvpyl0pX3b0KDz4IJw9axP7VUZLnzsHzz0HTz0FNWrYMowmdv+lyV0pX3X6tL0x6cABmDcPqlS54qn799u6+uefQ//+duPqf/8760JVWU/LMkr5ovh4u1FpdDTMmmUb06/gp5/s4v7cOVi40I6ZUf5Pk7tSviY52ba4LFlil+LNml3x1DVroGlTW46PjLzq4l75GS3LKOVLjLF3nE6ebIe+PPvsFU+dNw8aNrTTHNeu1cSe3WhyV8qXvP++3ZG6a1e7Hd4VhIfb/U0rV7ar92vsmqf8kFvJXUS6icg2EdkuIt1dx4qIyDIR2eN6LOyZUJXK5r74wo5nbN3aJvgrjBUYNszenFS/vh38pWMEsqcMJ3cRqQI8D9QEqgJNRKQC0A+IMMZUACJcz5VS7pg/3/YxNmxok3waYwWMsbm/d297rXXBAsifP+tDVd7BnZV7JWC9MeasMSYRWAU8AjQHwl3nhAMt3AtRqWxu3Tq7x121anYD05w5/3FKYqLd33ToUHuD0tSpdotUlX25k9y3AfVEpKiI5AEeBm4CShhjYgFcj2n+UigiHUQkSkSidLclpa5g+3a7G3Xp0raPMY2l+Llz8NhjdgDYG2/AmDEQEOBArMqrZLgV0hizU0TeA5YBp4GtQOJ1vH8cMA4gNDRUd25WKrWDB+1YgaAgWLrUtr2kcvy47YRcswY++gg6d3YgTuWV3Lqgaoz53BhT3RhTD/gT2AMcFpGSAK7HI+6HqVQ2c+yYnQ1w8qQdAJNGu0tsLNxzD6xfb+ewa2JXKbnbLVPc9fgf4FFgCjAXuLinVztgjjufoVS2c+aMHSvw66+2WT04+B+n7N0Ld99t9zidPx9atXIgTuXV3L1DdaaIFAUSgE7GmL9EZAgwXUTaA78BYe4GqVS2kZBgL55u2mQvnqaxz93mzbZak5RkZ8TUrOlAnMrruZXcjTF10zh2DLjfne+rVLaUnGzbHRcuhHHjoMU/G81WrrQ19kKFbBm+YsWsD1P5Br1DVSlv0a+f3T3jnXfg+ef/8fKsWXYAWOnStjtSE7u6Gk3uSnmD4cPtaIFOneCVV/7x8mefQVgYVK9uB4CVLu1AjMqnaHJXymmTJtlhYGFhMHLkZWMFjLHzwTp0sM0z330HRYs6GKvyGZrclXLSokV2suN999kkn+Luo+Rk6NHDLuTbtIG5c+3WeEqlhyZ3pZyyYYMdAhMcDLNnXzYvID4e2ra1C/lu3WwpPjDQwViVz9HNOpRyws6d8PDDULKk7Y4pUODSS2fO2Jy/eLEtyfTrd8UBkEpdkSZ3pbJaTIxtewkMtLsplShx6aVjx+z9S5s22Yuozz3nYJzKp2lyVyor/fmnTezHj8OqVXDLLZdeOnjQvvTLLzBjht1sQ6mM0uSuVFY5e9ZuaLp3r625VKt26aVdu2w3zPHj9qX69Z0LU/kHTe5KZYXERGjZEr7/HqZPh3vvvfTSxo22/B4QYBfzKXK+Uhmm3TJKZTZjbKP6/Pl22Pr//nfppaVLbRdkgQJ2E2tN7MpTNLkrldkGDPh7J42OHS8dnjrV7sNxyy02sZcv72CMyu9oclcqM334IQwZAi+8YJO7y+jR8MQTUKuWLcWULOlgjMovaXJXKrN8/bW9xfTRR+Hjj0EEY2yO79LFXltdssROeFTK0/SCqlKZYelSaNfOtr1MngwBASQl2aT+ySfwzDN2qm8O/RuoMomu3JXytI0b7Wq9cmX49lsICuLCBWjd2ib2vn3h8881savMpf95KeVJu3fbW0yLF7dDwQoW5NQpe0NSRAQMGwa9ejkdpMoONLkr5SmHDtlbTEVsMb1kSeLi4KGHYMsWCA+Hp55yOkiVXWhyV8oTjh+3G5seO2b3wqtQgf37ba4/eBDmzLELeqWyiiZ3pdx17pzd2HTXLjvhsUYNtm2zif3sWVi2DO6+2+kgVXajF1SVckdior1SumYNfPUVNGjA2rVQ17V1fGSkJnblDE3uSmWUMfDii7bmMnIkPP44CxZAw4ZQrJi967RKFaeDVNmVW8ldRHqIyHYR2SYiU0QkSESKiMgyEdnjeizsqWCV8iqvvQbjx9t98Lp04csvoXlzuP12u5AvW9bpAFV2luHkLiKlgK5AqDGmChAAtAL6ARHGmApAhOu5Uv7lo4/g3XftbhrvvMOIEX/fs7Rihe2EVMpJ7pZlcgC5RSQHkAc4BDQHwl2vhwMt3PwMpbzLtGl2Y9MWLTBjPqFff6FXLzvsccECyJ/f6QCVciO5G2N+B4YBvwGxwAljzFKghDEm1nVOLJDmGkZEOohIlIhExcXFZTQMpbLWd9/Znavr1CHxy695rmMO3nvPDnucOvWyPa6VcpQ7ZZnC2FV6OeDfQF4ReTK97zfGjDPGhBpjQosVK5bRMJTKOtHR9lbTihU5N20ujz2ZmwkT7CCwMWPsZhtKeQt3+twbAL8aY+IARGQW8F/gsIiUNMbEikhJ4IgH4lTKWXv22FtNixbl+LQlNGtZiDVrbOm9c2eng1Pqn9ypuf8G1BKRPCIiwP3ATmAu0M51TjtgjnshKuWw2Fh7R5IxxH4VwT2tSrJ+vZ3oq4ldeasMr9yNMRtEZAbwA5AIbAbGAfmA6SLSHvs/gDBPBKqUI06csCv2I0fYN2kdD7S7hcOH7Y55DzzgdHBKXZlb4weMMW8Ab6Q6fAG7ilfKt507ZxvXt29n88jVPPRiMImJsHw51KzpdHBKXZ3eoapUWi5csDPZV69mZd9F1O9fm5w57c1JmtiVL9DkrlRqCQl2XszixXzbcTGNhjWgVClYtw4qVnQ6OKXSR5O7UiklJdlbTWfP5qPHVvLYpw9QrZodAFa6tNPBKZV+mtyVuig5GZ5/nqQp0+h29ya6zryHpk3tfUtFizodnFLXR5O7UmAnPHbpwpmJ03j0tu2MWhtKjx4wcybkzet0cEpdP92sQyljoHdvDo2ZTdMSu9mypxSjR0OnTk4HplTGaXJX6s03+Wn4Ehrn28afpwszd67olnjK52lyV9nbkCEseXs9YYEbyV8giDULhJAQp4NSyn1ac1fZ16hRfNr/VxrLQm6+PYgNGzSxK/+hyV1lS8mffkafbufpyKc82EiIjBRtdVR+RZO7ynbOTZjC4x0L8z59eOmFJObMvUE32FB+R2vuKls5PH4ezZ8vx0ZqMuK9eLr3zomI01Ep5Xma3FW2sXPsKh5+sQqHbyjJrK8u0KJ1bqdDUirTaHJX2cLy96N5tE9VgnIksmppInfem8/pkJTKVFpzV37viwG7ebBPMKVzHWXDpgBN7Cpb0OSu/JYx8Nqzv/PM4Nuon2cTa7cVpExIYafDUipLaHJXfun8eXiy8Z8MnFiK9vmnsXB7GQqW143YVfahyV35naNHoeHdZ/l6UREGFxzCZ1vvIrBsKafDUipL6QVV5Vf27IGHG8Zz8MANTCvYgcc39YZyZZ0OS6kspyt35TciI6FWzSSOHzzF8oKP8vi67lChgtNhKeUITe7KL3z9NTRoYCh2Zj/r8zbgvyvehdtvdzospRyjyV35NGNg4EBo0wZqB2xiXc57uWXZWKhWzenQlHJUhpO7iNwmIltS/DkpIt1FpIiILBORPa5H7T1TmSI+Hp55Bl57DdoWmscS8wBFFn4Fd93ldGhKOS7Dyd0Ys9sYE2KMCQFqAGeB2UA/IMIYUwGIcD1XyqP++gsaNYLwcHjz358SfuZ/5Jr7DdSr53RoSnkFT5Vl7gf2GWMOAM2BcNfxcKCFhz5DKQB++QX++19Ys8YwqcLbvHGkMzJzBjRs6HRoSnkNT7VCtgKmuL4uYYyJBTDGxIpI8bTeICIdgA4A//nPfzwUhvJ369dDs2aQmGhYVqUn92wdBVOnQtOmToemlFdxe+UuIjmBZsA31/M+Y8w4Y0yoMSa0WDG9c1Bd24wZcO+9kD+f4fsqHbhny0hblwkLczo0pbyOJ8oyDwE/GGMOu54fFpGSAK7HIx74DJWNGQPvv29zePVqhvWVnuG2yPEwdiw8+aTT4SnllTyR3Fvzd0kGYC7QzvV1O2COBz5DZVMJCdCxI/TpAy0fTybiP89QbGE4fPghdOjgdHhKeS23kruI5AEaArNSHB4CNBSRPa7XhrjzGSr7OnkSmjSBceNgQH/D1/leIGhaOAweDN26OR2eUl7NrQuqxpizQNFUx45hu2eUyrDffrOJfedOGP+Zof3WrjBhvG1q76fdtUpdi96hqrxOdLS9D+nAAVi00ND+574wejT06gVvveV0eEr5BE3uyqvMnWvvQ8qVC9atgwbr3rZXU1980T7qbtZKpYsmd+U1Ro2CFi2gcmXbz155wVB48014+mm7ctfErlS6aXJXjktKgq5d7TXSFi1g5Ur41zcfQd++0LIljB8PN+h/qkpdD/0boxx1+rRN6B99ZEvq33wDeb4eb7N98+YwaRIEBDgdplI+R3diUo45dMh2xGzdCmPG2LI6kyfb/vVGjWDaNAgMdDpMpXySJnfliK1bbWI/fhzmzYOHHwZmzoR27aB+fZg1y15VVUpliJZlVJZbvBjq1LFjBdascSX2BQugdWuoWdO2zOTO7XSYSvk0Te4qS40da1fs5cvDhg1QtSoQEQGPPQbBwbBoEeTL53SYSvk8Te4qSyQnw8sv27p6o0awejWUKoVdujdrZjeyXrIEChZ0OlSl/ILW3FWmO3sW2ra1ZfROnezMrxw5gE2bbE2mdGn47jsoWvSa30splT6a3FWmOnzYLsw3bYIPPrC97CLYK6oPPgg33mjLMiVKOB2qUn5Fk7vKNDt22IV5XBzMnm3b1gE7DaxhQ8ib1yb20qUdjVMpf6Q1d5UpIiLsPqcXLsCqVSkS+969cP/99o7TiAgoV87ROJXyV5rclUcZA599Zi+a3nSTnRETGup68cABm9jj422N/dZbHY1VKX+myV15zKlT9sJphw5w3322EaZMGdeLhw7ZxH7iBCxdClWqOBqrUv5Ok7vyiC1boEYNmDIF3nkHFi5M0dV45IhN7IcP2zuYqld3NFalsgO9oKrcYgx88gn07Gk7GVessPPYL/nzT3jgAdfOG4ugVi3HYlUqO9GVu8qw48chLMz2rt93n129X5bYT560xfedO+Hbb+GeexyLVansRpO7ypCNG6FaNZgzB4YOhfnzoVixFCecOQONG8PmzTBjhl29K6WyjCZ3dV2MgREj4O677UiB1auhd+9Ue2mcP297H9etsyN8mzZ1LF6lsitN7irdjh2zd5v26mWHf23ZArVrpzrpjz/snUsRETBxIjz+uCOxKpXduZXcRaSQiMwQkV0islNEaotIERFZJiJ7XI+FPRWscs6aNRASYrsYR42yc2IKp/43u3SpHfO4fj18+SU89ZQjsSql3F+5jwQWG2MqAlWBnUA/IMIYUwGIcD1XPio5GQYNsvtn5MplKy1duqTaqzohAfr1s7NiihWzg2TatnUqZKUUbiR3ESkA1AM+BzDGxBtjjgPNgXDXaeFAC3eDVM44fNg2u7zyiu2K+eEH28t+mV9/tS0y770HL7xgr7RWruxIvEqpv7mzcr8ZiAMmishmERkvInmBEsaYWADXY/G03iwiHUQkSkSi4uLi3AhDZYaICFthiYy04wS+/hoKFEh10jff2JaZHTtg+nS7E0eePI7Eq5S6nDvJPQdQHfjEGFMNOMN1lGCMMeOMMaHGmNBil/XQKSclJsLrr9uhjUWK2IX4c8+lKsOcOwcdO9qLpRUr2iurYWGOxayU+id3knsMEGOM2eB6PgOb7A+LSEkA1+MR90JUWeX33+2UgHfegaeftqXzO+5IddL27Xaf008/hb597dJeJzsq5XUynNyNMX8AB0XkNteh+4EdwFygnetYO2COWxGqLLFwoS3DREfbRpcJE+y49UuMgfHj4c477ayYJUtgyBAIDHQsZqXUlbk7W6YLMFlEcgK/AM9g/4cxXUTaA78B+vu6F0tIgAEDYNgwuz/19Olw222pTjpxwl4snTYNGjSASZPgX/9yJF6lVPq4ldyNMVuA0DReut+d76uyxv790KoVbNhgN64ePhxy50510saN9qTffoPBg6FPn1S3oyqlvJH+Lc2mZs2yjS47d9rV+pgxqRJ7crJdzl+cMxAZaXvZNbEr5RP0b2o2c/68vcZtChYAAA2cSURBVAnpscegfHk71+sfjS5HjtihX7172xkxmzenMWdAKeXNNLlnI3v22H1NR4+GHj1g7Vq4+eZUJ11scF+xwg5q/+abNOYMKKW8nSb3bGLKFLsB0v79dkzviBGQM2eKExIT4dVXbYN74cK21t6xY6oGd6WUr9Dk7ufOnrV7mj7xhO2G2bLFTna8zG+/2eEx774Lzz5rG9yDg50IVynlIbrNnh/bsQNatoRt26B/f3jrrTTa0mfPtgk9KcnOGGjd2pFYlVKepSt3P2QMfPGFvd/o4p7UgwalSuznz9v98R599O8rq5rYlfIbmtz9zOnT0K4dPPMM3HUXbN1qJ/FeZtcu++KYMXbnjbVr4ZZbHIlXKZU5NLn7ka1b7UjeyZNtCWbZMihZMsUJF5f0NWrAoUOwYIHtZb/syqpSyh9ocvcDxthpu3fdBadO2W7G11+HgIAUJ506ZTfQSLmkf/hhx2JWSmUuTe4+7sQJe9H0xRfh3nttN0z9+qlOio62fZBTptiRj8uWwb//7US4Sqksosndh23aZEcIzJplN0JasACKp9waxRj48EN7d+n587Bype1lv2xJr5TyR5rcfZAx8MEHduxLUhKsXp3GPK+jR21De48etvyydSvUretYzEqprKV97j7m2DFbNp83z459mTDB7ph0mZUroU0bm+A/+si2POqdpkplK7py9yFr10JIiO1bHznS3n90WWJPTIQ33oD77oN8+ews386dNbErlQ1pcvcBycl2lPo999iuxXXroGvXVDk7Jsbukff22/DUU/YiakiIYzErpZylZRkvd/iwzdVLl9r9qMeNg4IFU500b57d9PTCBbtHXtu2ToSqlPIiunL3YsuX28X36tV2P+qpU1Ml9gsXoHt3e+G0TBn44QdN7EopQJO7V0pKsqXzBg2gUCFbOu/QIVUZ5uefbYvjyJHQrRt8/z3ceqtjMSulvIuWZbzM77/bRpdVq+yMmNGj7bXRy0yaZO9aypXLDmf/xwxfpVR2pyt3L7JokS3DbNpkR8B88UWqxH5xKthTT9n5MFu3amJXSqVJk7sXOH7ctqI//LAd9BUdbXP4ZbZssQn9q69szWb5cihd2pF4lVLez63kLiL7ReQnEdkiIlGuY0VEZJmI7HE96gacV2CMzdW33WYHf3XtauvrFSumOmn0aDvs6/RpOxXszTd1hIBS6qo8sXK/1xgTYowJdT3vB0QYYyoAEa7nKpVdu2xbetu2ULasLcWMHAm5c6c46dgxeOQR6NLF7m26dWsaU8GUUuqfMqMs0xwId30dDrTIhM/wWWfPwiuv2C1KN2+GTz6xNyVVr57qxMhIW4BfuNAOkpk3D2680ZGYlVK+x93kboClIhItIh1cx0oYY2IBXI/F03qjiHQQkSgRiYqLi3MzDN+wYAFUrmy3vGvVyq7eO3ZMVWFJSrJjeevXh6Ag2+LYvbuOEFBKXRd3WyHvNsYcEpHiwDIR2ZXeNxpjxgHjAEJDQ42bcXi1gwdtK/rs2VCpEqxYcYXqyqFDtg/y4uCvTz6B/PmzOFqllD9wa+VujDnkejwCzAZqAodFpCSA6/GIu0H6qoQEu4tdpUp/b1Kd5mYaCQkwcSJUrQobN9qvJ03SxK6UyrAMJ3cRySsi+S9+DTwAbAPmAhcb+doBc9wN0hetXWvr6L172yGNO3ZA//6ptiu9cMEOi7n1Vnj2WXtlNTrazonRMoxSyg3urNxLAGtEZCuwEVhgjFkMDAEaisgeoKHrebZx9KjN03Xq2C3wvv0W5s61efuS8+fh44+hfHl44QW7fdK8eXbVflkfpFJKZUyGa+7GmF+AqmkcPwbc705Qvig52W6c0bcvnDxpd0Z6/XXImzfFSWfP2pX60KEQG2u3Uvr8c9vmqCt1pZQH6WwZD/jxR9v18v33die7MWOgSpUUJ5w+bQ8OHw5HjtidrCdPtsV3TepKqUyg4wfccOoU9Opla+t79thZMKtWpUjsJ07Au+/acbx9+9q+9chIOzrg3ns1sSulMo2u3DPAGJg507af//67Hcc7eHCKLe/+/NPebjpqlB0c06QJvPqqHSGglFJZQJP7ddq3z25Lunix7VycMQNq1XK9ePQojBhhZ8GcOmVHB7z6ahq3nyqlVObS5J5OFy7Y66CDBkGOHHYiQOfO9mv++MPW08eMgXPnICzMJvU77nA6bKVUNqXJPR0iIuCll+zmR2FhNrGXKoWtyQwdajtg4uPhiSdgwAB715JSSjlIL6heRWyszdcNGtiRL4sWwfTpUCrxgM32N99sV+utW9tBMZMmaWJXSnkFXbmnISnJ5uxXX7X3G73+OvTrB7kP7YPnBkN4uO10eeYZ+0K5ck6HrJRSl9HknsqmTbZn/Ycf7L1FH38MFZJ3Q8dBtjc9Rw57Qp8+cNNNToerlFJp0uTucvy4LZePHQv/+hdMnQqPV96OvPEuTJtmN6Pu2tUOiylZ0ulwlVLqqrJ9cjfGLsh79bKdjF26wNv/+5GCI9+GVjPt/IDevaFnTzsDRimlfEC2Tu47d9rroitXQs2asGj4Dqp/0x/qzYUCBWzRvXt3KFrU6VCVUuq6ZMvkfvYsDBxoZ63nzQtj++zjuR+7EdB2ARQuDG+9ZUswhQo5HapSSmVItkvu8+fb0sv+/fDUg3/w/tnOFB860+5POniwXcoXKOB0mEop5ZZs0+f+2292GkDTppA7+TQrg7sSvqQkxX9eY5fw+/fbtkZN7EopP+D3K/eEBPjwQ3jzTYNJSmbwTZ/S87fu5Ewqbgd7Pfcc5M7tdJhKKeVRfp3c16yBF180bNsmNC0YyaizT1FWDHwyyt6AlCuX0yEqpVSm8MvkHhcHffsYJn4h/Ccwlm95keZFt8Hw16Bt21QbmSqllP/xq+SenAyff5ZMv5cTOHn6BvoynNfKTCXv672g9UzXCEellPJ/fpPttkYn8mKrP/l+b3HqsZ4xt4yg8sDWEBYNAQFOh6eUUlnK55P7qWPxvNFqN6O+q0RhhC9ueo2nRoQgj86GG7JNM5BSSl3Gp5N7VPh2WrQvwu9Jd9Ch2CwGf5iHIq3f1r1JlVLZntvJXUQCgCjgd2NMExEpAkwDygL7gceNMX+5+zlpuaVeKW4vtI8Zr+ynVvdHNKkrpZSLJ+oW3YCdKZ73AyKMMRWACNfzTFG4XCGWHq1BrR61NbErpVQKbiV3ESkNNAbGpzjcHAh3fR0OtHDnM5RSSl0/d1fuHwJ9gOQUx0oYY2IBXI9pzskVkQ4iEiUiUXFxcW6GoZRSKqUMJ3cRaQIcMcZEZ+T9xphxxphQY0xosWLFMhqGUkqpNLhzQfVuoJmIPAwEAQVE5CvgsIiUNMbEikhJ4IgnAlVKKZV+GV65G2P6G2NKG2PKAq2A5caYJ4G5QDvXae2AOW5HqZRS6rpkxl0+Q4CGIrIHaOh6rpRSKgt55CYmY8xKYKXr62PA/Z74vkoppTJG789XSik/JMYYp2NAROKAA258ixuBox4Kx0n+8nOA/izeyF9+DtCf5aIyxpg02w29Irm7S0SijDGhTsfhLn/5OUB/Fm/kLz8H6M+SHlqWUUopP6TJXSml/JC/JPdxTgfgIf7yc4D+LN7IX34O0J/lmvyi5q6UUupy/rJyV0oplYImd6WU8kM+ndxFpJGI7BaRvSKSaZuCZDYRmSAiR0Rkm9OxuEtEbhKRFSKyU0S2i0g3p2PKCBEJEpGNIrLV9XO85XRM7hKRABHZLCLznY7FHSKyX0R+EpEtIhLldDwZJSKFRGSGiOxy/X2p7dHv76s1d9f2fj9j59fEAJuA1saYHY4GlgEiUg84DXxpjKnidDzucE0CLWmM+UFE8gPRQAtf+/ciIgLkNcacFpFAYA3QzRiz3uHQMkxEegKhQAFjTBOn48koEdkPhBpjfPomJhEJByKNMeNFJCeQxxhz3FPf35dX7jWBvcaYX4wx8cBU7C5QPscYsxr40+k4PMEYE2uM+cH19SnsFoylnI3q+hnrtOtpoOuPb66EuOKuacohIlIAqAd8DmCMifdkYgffTu6lgIMpnsfgg0nEn4lIWaAasMHZSDLGVcbYgt2TYJkxxid/Dpe0dk3zVQZYKiLRItLB6WAy6GYgDpjoKpWNF5G8nvwAX07uae2I7bMrK38jIvmAmUB3Y8xJp+PJCGNMkjEmBCgN1BQRnyyZubtrmhe62xhTHXgI6OQqa/qaHEB14BNjTDXgDODR64a+nNxjgJtSPC8NHHIoFpWCq0Y9E5hsjJnldDzucv26vBJo5HAoGXVx17T92PLlfa5d03ySMeaQ6/EIMBtbovU1MUBMit8GZ2CTvcf4cnLfBFQQkXKuixGtsLtAKQe5LkR+Duw0xoxwOp6MEpFiIlLI9XVuoAGwy9moMuYqu6b5HBHJ67pQj6uM8QDgc11mxpg/gIMicpvr0P2AR5sOPLJZhxOMMYki0hlYAgQAE4wx2x0OK0NEZApQH7hRRGKAN4wxnzsbVYbdDbQFfnLVqwEGGGMWOhhTRpQEwl1dWTcA040xPt1C6CdKALPtGoIcwNfGmMXOhpRhXYDJrsXpL8AznvzmPtsKqZRS6sp8uSyjlFLqCjS5K6WUH9LkrpRSfkiTu1JK+SFN7kop5Yc0uSullB/S5K6UUn7o/1n5VVY8onGqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# GM(2, 1)模型 (DGM(2, 1)模型是它的一个简化, 实现相似, 这里不再实现)\n",
    "# 例15.3\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import sympy\n",
    "\n",
    "# 1.级比检验\n",
    "x0 = np.array([41, 49, 61, 78, 96, 104], dtype=np.float)\n",
    "lamb = np.array([x0[i] / x0[i+1] for i in range(x0.shape[0]-1)])\n",
    "print('级比:', lamb)\n",
    "Theta = (np.exp(-2/(x0.shape[0]+1)), np.exp(2/(x0.shape[0]+2)))\n",
    "print('可容覆盖:', Theta)\n",
    "OK = [i <= Theta[1] and i >= -Theta[0] for i in lamb]\n",
    "print('可容性:', OK)\n",
    "\n",
    "# GM(2, 1)建模\n",
    "x1 = np.array([x0[:i].sum() for i in range(1, x0.shape[0]+1)])\n",
    "print('1-AGO序列:', x1)\n",
    "alpha1_x0 = np.array([x0[i] - x0[i-1] for i in range(1, x0.shape[0])])\n",
    "print('1-IAGO序列:', alpha1_x0)\n",
    "z1 = np.array([(x1[i]+x1[i+1])/2 for i in range(x1.shape[0]-1)])\n",
    "print('均值生成序列:', z1)\n",
    "B = np.array(np.c_[-x0[1:], -z1, np.ones(z1.shape[0])])\n",
    "print('B矩阵:', B, sep='\\n')\n",
    "Y = alpha1_x0\n",
    "print('Y值:', Y)\n",
    "u = np.matmul(np.linalg.inv(np.matmul(B.T, B)), np.matmul(B.T, Y))\n",
    "print('最小二乘法得到a1, a2, b:', u)\n",
    "\n",
    "sympy.init_printing()\n",
    "# 定义符号常量x 与 f(x)\n",
    "t = sympy.Symbol('t')\n",
    "a1, a2, b = sympy.symbols('a1, a2, b')\n",
    "x1_t = sympy.Function('x1_t')\n",
    "# 微分方程ode\n",
    "ode = sympy.Eq(x1_t(t).diff(t, t) + a1*x1_t(t).diff(t) + a2*x1_t(t), b)\n",
    "init1 = sympy.Eq(x1_t(0), x1[0])\n",
    "init2 = sympy.Eq(x1_t(1), x1[1])\n",
    "# 调用dsolve函数,返回一个Eq对象，hint控制精度\n",
    "ode_sol = sympy.dsolve((ode), x1_t(t))\n",
    "print('得到带参符号解:', ode_sol)\n",
    "# 初始条件字典\n",
    "x1_to_value = {x1_t(0): 41, x1_t(x1.shape[0]-1): x1[-1]}\n",
    "# 确定待定常数集\n",
    "known_params = {a1, a2, b, t}\n",
    "free_params = ode_sol.free_symbols - set(known_params)\n",
    "# 将初始条件带入微分方程的后得到的多元方程\n",
    "C_eqs = [sympy.Eq(ode_sol.lhs.subs(t, 0).subs(x1_to_value), ode_sol.rhs.subs(t, 0)),\n",
    "         sympy.Eq(ode_sol.lhs.subs(t, x1.shape[0]-1).subs(x1_to_value), ode_sol.rhs.subs(t, x1.shape[0]-1))\n",
    "        ]\n",
    "print('待定常数: ', free_params)\n",
    "print('关于待定常数的多元方程组:', C_eqs, sep='\\n')\n",
    "sol_params = sympy.solve(C_eqs, free_params)\n",
    "print('解得待定常数:', sol_params, '\\n')\n",
    "known_params_value = {a1: u[0], a2: u[1], b: u[2]}\n",
    "final_sol_eq = sympy.Eq(ode_sol.lhs, ode_sol.rhs.subs(sol_params).subs(known_params_value))\n",
    "print('最终解:', final_sol_eq.lhs, '=', final_sol_eq, '\\n')\n",
    "\n",
    "# 预测\n",
    "x1_pred = np.array([final_sol_eq.rhs.subs(t,i) for i in range(x1.shape[0])])\n",
    "print('x1预测值:', x1_pred)\n",
    "x0_pred = np.array([x1_pred[i] - x1_pred[i-1] if i != 0 else x1_pred[i] for i in range(x1_pred.shape[0])])\n",
    "print('x0预测值:', x0_pred)\n",
    "\n",
    "# 检验\n",
    "test_table_columns = ['Raw Data', 'Predction', 'Error', 'Relative_Error']\n",
    "test_table = pd.DataFrame(np.c_[x0,\n",
    "                                x0_pred,\n",
    "                                x0 - x0_pred,\n",
    "                                np.divide((x0 - x0_pred), x0)],\n",
    "                          columns=test_table_columns)\n",
    "test_table['Predction'] = test_table['Predction'].apply(lambda x: format(float(x), '.4f'))\n",
    "test_table['Error'] = test_table['Error'].apply(lambda x: format(float(x), '.4f'))\n",
    "test_table['Relative_Error'] = test_table['Relative_Error'].apply(lambda x: format(float(x), '.2%'))\n",
    "print('\\n', test_table, sep='')\n",
    "\n",
    "plt.plot(np.linspace(0, x0.shape[0]-1, x0.shape[0]), x0, c='red', label='Raw Data')\n",
    "x1_pred = np.array([final_sol_eq.rhs.subs(t,i) for i in range(7)])\n",
    "x0_pred = np.array([x1_pred[i] - x1_pred[i-1] if i != 0 else x1_pred[i] for i in range(x1_pred.shape[0])])\n",
    "plt.plot(np.linspace(0, x0_pred.shape[0]-1, x0_pred.shape[0]), x0_pred, c='blue', label='Pred Data')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "一次累加序列: [ 4.93  7.26 11.13 15.48 22.11 29.26 34.63 41.02 48.83 57.18]\n",
      "均值生成序列: [ 6.095  9.195 13.305 18.795 25.685 31.945 37.825 44.925 53.005]\n",
      "B矩阵:\n",
      "[[  -6.095      37.149025]\n",
      " [  -9.195      84.548025]\n",
      " [ -13.305     177.023025]\n",
      " [ -18.795     353.252025]\n",
      " [ -25.685     659.719225]\n",
      " [ -31.945    1020.483025]\n",
      " [ -37.825    1430.730625]\n",
      " [ -44.925    2018.255625]\n",
      " [ -53.005    2809.530025]]\n",
      "Y值: [2.33 3.87 4.35 6.63 7.15 5.37 6.39 7.81 8.35]\n",
      "最小二乘法得到a, b: [-0.35762668 -0.00410372]\n",
      "x1预测值: [  4.93         7.04462985  10.06840608  14.39219957  20.57492891\n",
      "  29.41580956  42.05766501  60.13465479  85.98351544 122.94562055]\n",
      "x0预测值: [ 4.93        2.11462985  3.02377624  4.32379348  6.18272934  8.84088065\n",
      " 12.64185545 18.07698979 25.84886065 36.96210511]\n",
      "\n",
      "   Raw Data Predction     Error Relative_Error\n",
      "0      4.93    4.9300    0.0000          0.00%\n",
      "1      2.33    2.1146    0.2154          9.24%\n",
      "2      3.87    3.0238    0.8462         21.87%\n",
      "3      4.35    4.3238    0.0262          0.60%\n",
      "4      6.63    6.1827    0.4473          6.75%\n",
      "5      7.15    8.8409   -1.6909        -23.65%\n",
      "6      5.37   12.6419   -7.2719       -135.42%\n",
      "7      6.39   18.0770  -11.6870       -182.89%\n",
      "8      7.81   25.8489  -18.0389       -230.97%\n",
      "9      8.35   36.9621  -28.6121       -342.66%\n"
     ]
    }
   ],
   "source": [
    "# 灰色Verhulst预测模型\n",
    "# 例15.5\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# 原始数据\n",
    "x0 = np.array([4.93, 2.33, 3.87, 4.35, 6.63, 7.15, 5.37, 6.39, 7.81, 8.35], dtype=np.float)\n",
    "\n",
    "# 计算1次累加序列\n",
    "x1 = np.array([x0[:i].sum() for i in range(1, x0.shape[0]+1)])\n",
    "print('一次累加序列:', x1)\n",
    "z1 = np.array([(x1[i]+x1[i+1])/2 for i in range(x1.shape[0]-1)])\n",
    "print('均值生成序列:', z1)\n",
    "B = np.array([[-z1[i], z1[i]*z1[i]] for i in range(x1.shape[0]-1)])\n",
    "print('B矩阵:', B, sep='\\n')\n",
    "Y = np.array([x0[i] for i in range(1, x0.shape[0])])\n",
    "print('Y值:', Y)\n",
    "u = np.matmul(np.linalg.inv(np.matmul(B.T, B)), np.matmul(B.T, Y))\n",
    "print('最小二乘法得到a, b:', u)\n",
    "a, b = u[0], u[1]\n",
    "x1_k = lambda k: (x0[0] - b/a) * np.exp(-a*k) + b/a\n",
    "\n",
    "# 预测\n",
    "x1_pred = np.array([x1_k(k) for k in range(x0.shape[0])])\n",
    "print('x1预测值:', x1_pred)\n",
    "x0_pred = np.array([x1_pred[i] - x1_pred[i-1] if i != 0 else x1_pred[0] for i in range(x1_pred.shape[0])])\n",
    "print('x0预测值:', x0_pred)\n",
    "\n",
    "# 检验\n",
    "test_table_columns = ['Raw Data', 'Predction', 'Error', 'Relative_Error']\n",
    "test_table = pd.DataFrame(np.c_[x0,\n",
    "                                x0_pred,\n",
    "                                x0 - x0_pred,\n",
    "                                np.divide((x0 - x0_pred), x0)],\n",
    "                          columns=test_table_columns)\n",
    "test_table['Predction'] = test_table['Predction'].apply(lambda x: format(float(x), '.4f'))\n",
    "test_table['Error'] = test_table['Error'].apply(lambda x: format(float(x), '.4f'))\n",
    "test_table['Relative_Error'] = test_table['Relative_Error'].apply(lambda x: format(float(x), '.2%'))\n",
    "print('\\n', test_table, sep='')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
